rm(list=ls())
library(statnet)
#setwd("C:\\Documents and Settings\\Laura Metzger\\Desktop\\Dissertation\\Chapter2\\Data\\")
data<-read.csv("NH07.csv")
load("NHnet07.RData")
diag(AM)<-0

AMnet<-network(AM)
r<-as.vector(AM)
r<-r[r!=0]
set.edge.attribute(AMnet, "values", r, e=1:length(r))
AMnet<-network.adjacency(AM, AMnet, ignore.eval=FALSE, names.eval="values")
get.edge.attribute(AMnet$mel, "values")

same.party<-matrix(0, nrow(AM), nrow(AM))
for(i in 1:nrow(AM)){
for(j in 1:nrow(AM)){
same.party[i,j]<-ifelse(data$party[i]==data$party[j], 1, 0)
}
}


same.dist<-matrix(0, nrow(AM), nrow(AM))
for(i in 1:nrow(AM)){
for(j in 1:nrow(AM)){
same.dist[i,j]<-ifelse(data$dist[i]==data$dist[j], 1, 0)
}
}
same.dist<-ifelse(is.na(same.dist)==TRUE, 0, same.dist)

diag(same.dist)<-0

vec1<-as.vector(same.party)
vec2<-as.vector(same.dist)
vec3<-vec1*vec2
interact<-matrix(vec3, ncol=nrow(AM))

mmd<-vector(mode="numeric", length=nrow(AM))
for(i in 1:nrow(same.dist)){
mmd[i]<-ifelse(sum(same.dist[i,])>1 & sum(same.party[i,which(same.dist[i,]==1)]==0)>0, 1, 0)
}
set.vertex.attribute(AMnet, "mmd", mmd, v=1:nrow(AM))

mod<-ergmm(AMnet~latentcov(same.party)+latentcov(same.dist)+latentcov(interact)+receivercov("mmd")+latent(d=2), family="Poisson", response="values", verbose=1, control=ergmm.control(sample=10000, burnin=10000))

dput(mod, file="resultsNH.txt")
summary(mod)

r<-predict.ergmm(mod)
f<-AM-r

pdf("DensityDifferenceNH.pdf")
plot(density(f), main="Density of Differences between Observed and Predicted Connections", xlab="Difference in Connections")
dev.off()

pdf("DensityObsPreNH.pdf")
plot(density(AM), col="red", xlim=c(-5, 20), main="Density of Observed and Predicted Connections in the MD House")
lines(density(r), col="blue", xlim=c(-5, 20))
legend(10, 0.25, c("Observed","Predicted"), col=c("red", "blue"), lty=c(1,1))
dev.off()

pdf("positionsNH.pdf")
plot(mod, vertex.col=data$party+1, main="MKL Latent Positions of MD House in 2007", label=TRUE, edge.col=0, what="mkl")
abline(0,0)
abline(v=0)
legend(-5, 6, c("Democrats","Republicans"), col=c("red","black"), lty=c(1,1), lwd=c(2,2))
dev.off()
